#' ---
#' title: "Analysis of Conjoint Experiment"
#' author: "Lukas F. Stoetzer"
#' date: "Feb 2022"
#' ---

  # Load Data
  df <- readRDS("dat_cnj.RDS")

# Model AMCE ==========

  # Estimation Function AMCE
  est_lm <- function(df){
    m <- lm(y ~ age + arrival + chance + child + job + gender + partisanship1 + partisanship2, data=df,weights = weight)
    clse_models <-  vcovCL(m, ~ ResponseId)
    d <- data.frame("est" = coef(m),
                    "se" = sqrt(diag(clse_models)),
                    "names" = names(coef(m))) %>%
               filter(names!="(Intercept)")
    return(d)
  } 
  
  
  # Estimation countrywise
  df_res <- df %>% 
    group_by(Cntry) %>%
    do(est_lm(.))
  
  # List of ACMEs
  df_names <- data.frame(
    names = c("age27","age42","age61", "age76",
              "jobNurse","jobPhysician","jobProfessor","jobUnemployed","jobCook",                
              "arrivalFirst","arrivalSecond","arrivalSame time",
              "chance20%","chance50%","chance80%",
              "partisanship1","partisanship2","partisanship3",
              "genderFemale","genderMale",
              "childYes","childNo"),
    attr   = c("Age","Age","Age", "Age",
               "Job","Job","Job","Job","Job",                  
               "Arrival","Arrival","Arrival",
               "Chance of Survival","Chance of Survival", "Chance of Survival",
               "Partisanship","Partisanship","Partisanship",
               "Gender","Gender",
               "Children","Children"),
    levl = c("27","42","61", "76",
             "Nurse","Physician","Professor","Unemployed","Cook",                  
             "First","Second","Same time",
             "20%","50%","80%",
             "Right","Left","None",
             "Female","Male",
             "Yes","No"))
  
  df_names <- crossing(df_names, data.frame("Cntry"=unique(df_res$Cntry)))

    # Create Baseline
  df_plot <- left_join(df_names, df_res) %>%
    mutate(est = ifelse(is.na(est),0,est),
           se = ifelse(is.na(se),0,se),
           baseline = ifelse(est == 0 & se ==0,"Baseline","ACME"),
           levl = as.factor(levl), 
           Cntry = factor(Cntry, levels = c("Countries combined","USA","Germany","Italy","Poland","Brazil"))) 
  
  df_plot$label <- as.factor(interaction(df_plot$levl, df_plot$attr))
  df_plot$position <- case_when(df_plot$label == "27.Age" ~ 1,
                                df_plot$label == "42.Age" ~ 2,
                                df_plot$label == "61.Age" ~ 3,
                                df_plot$label == "76.Age" ~ 4,
                                
                                df_plot$label == "First.Arrival" ~ 6,
                                df_plot$label == "Same time.Arrival" ~ 7,
                                df_plot$label == "Second.Arrival" ~ 8,
                                
                                df_plot$label == "20%.Chance of Survival" ~ 10,
                                df_plot$label == "50%.Chance of Survival" ~ 11,
                                df_plot$label == "80%.Chance of Survival" ~ 12,
                                
                                df_plot$label == "No.Children" ~ 14,
                                df_plot$label == "Yes.Children" ~ 15,
                                
                                df_plot$label == "Female.Gender" ~ 17,
                                df_plot$label == "Male.Gender" ~ 18,
                                
                                df_plot$label == "Cook.Job" ~ 20,
                                df_plot$label == "Nurse.Job" ~ 21,
                                df_plot$label == "Physician.Job" ~ 22,
                                df_plot$label == "Professor.Job" ~ 23,
                                df_plot$label == "Unemployed.Job" ~ 24,
                                
                                df_plot$label == "None.Partisanship" ~ 26,
                                df_plot$label == "Right.Partisanship" ~ 27,
                                df_plot$label == "Left.Partisanship" ~ 28)
  

  xaxis_ticks <- unique(df_plot$position)
  xaxis_lab <- gsub("\\."," ", unique(df_plot$label))
  
  vertival_lines <- c(5,9,13,16,19,25)
  

  # Figure
  ggplot(df_plot,aes(x=position,y=est,ymin=est-1.96*se,ymax=est+1.96*se, 
                     col=baseline)) + 
    geom_hline(yintercept = 0,alpha=0.3,col="red") +
    facet_wrap(~ Cntry) + theme_bw() +
    geom_pointrange(size=0.7) +
    scale_x_continuous(breaks=xaxis_ticks,
                       labels = xaxis_lab) +
    geom_vline(xintercept=vertival_lines, col="grey", linetype="dashed", size=1) +
    coord_flip() +  xlab("") + ylab("AMCE") + 
    scale_color_grey() +
    theme(legend.position = "none", text = element_text(size=12),
          #panel.grid.major.y =  element_blank(),
          panel.grid.minor.y =  element_blank())

  ggsave("appendix_figure11.pdf", height=9, width=12)

  
  

# Model Marginal Means ==========
  
  # Estimation Function MM
  est_mm <- function(df){

    # Create Partisianship
    df$partisanship_lab <- ifelse(df$partisanship1 == 1, "Right Partisianship",
                           ifelse(df$partisanship2 == 1, "Left Partisianship",
                                                         "None Partisianship"))
    
    # Linear Model
    m <- lm(y ~ age + arrival + chance + child + job + gender + partisanship_lab, 
            data=df,weights = weight)
    
    
    # Caclualte Marginal Means
    vars <- c("age", "arrival", "chance" ,"child" , "job" , "gender", "partisanship_lab")
    df_res <- do.call( rbind, lapply(vars, function(v) {  
                              df_mm <-  data.frame(emmeans(m,v)) 
                              names(df_mm)[1] <- "var"  
                              return(df_mm)
                            }))
    
    # Return
    return(df_res)
  } 
    
  # Estimate for all Countries
  df_res_mm <- df %>% 
    group_by(Cntry) %>%
    do(est_mm(.))

  # Prepare LAbels for Figure
  # List of ACMEs
  df_names <- data.frame(
    var = c("27","42","61", "76",
              "Nurse","Physician","Professor","Unemployed","Cook",                
              "First","Second","Same time",
              "20%","50%","80%",
              "Left Partisianship","Right Partisianship","None Partisianship",
              "Female","Male",
              "Yes","No"),
    attr   = c("Age","Age","Age", "Age",
               "Job","Job","Job","Job","Job",                  
               "Arrival","Arrival","Arrival",
               "Chance of Survival","Chance of Survival", "Chance of Survival",
               "Partisanship","Partisanship","Partisanship",
               "Gender","Gender",
               "Children","Children"),
    levl = c("27","42","61", "76",
             "Nurse","Physician","Professor","Unemployed","Cook",                  
             "First","Second","Same time",
             "20%","50%","80%",
             "Left","Right","None",
             "Female","Male",
             "Yes","No"))
  
  df_names <- crossing(df_names, 
                       data.frame("Cntry"=unique(df_res$Cntry)))
  
  
  # Create Plot Dataset
  df_plot <- left_join(df_names, df_res_mm) %>%
    mutate(levl = as.factor(levl), 
           Cntry = factor(Cntry, levels = c("Countries combined","USA","Germany","Italy","Poland","Brazil"))) 
  
  df_plot$label <- as.factor(interaction(df_plot$levl, df_plot$attr))
  df_plot$position <- case_when(df_plot$label == "27.Age" ~ 1,
                                df_plot$label == "42.Age" ~ 2,
                                df_plot$label == "61.Age" ~ 3,
                                df_plot$label == "76.Age" ~ 4,
                                
                                df_plot$label == "First.Arrival" ~ 6,
                                df_plot$label == "Same time.Arrival" ~ 7,
                                df_plot$label == "Second.Arrival" ~ 8,
                                
                                df_plot$label == "20%.Chance of Survival" ~ 10,
                                df_plot$label == "50%.Chance of Survival" ~ 11,
                                df_plot$label == "80%.Chance of Survival" ~ 12,
                                
                                df_plot$label == "No.Children" ~ 14,
                                df_plot$label == "Yes.Children" ~ 15,
                                
                                df_plot$label == "Female.Gender" ~ 17,
                                df_plot$label == "Male.Gender" ~ 18,
                                
                                df_plot$label == "Cook.Job" ~ 20,
                                df_plot$label == "Nurse.Job" ~ 21,
                                df_plot$label == "Physician.Job" ~ 22,
                                df_plot$label == "Professor.Job" ~ 23,
                                df_plot$label == "Unemployed.Job" ~ 24,
                                
                                df_plot$label == "None.Partisanship" ~ 26,
                                df_plot$label == "Right.Partisanship" ~ 27,
                                df_plot$label == "Left.Partisanship" ~ 28
                                )
  
  
  xaxis_ticks <- unique(df_plot$position)
  xaxis_lab <- gsub("\\."," ", unique(df_plot$label))
  
  vertival_lines <- c(5,9,13,16,19,25)
  
  # Figure
  ggplot(df_plot,aes(x=position,y=emmean,ymin=lower.CL,ymax=upper.CL)) + 
    geom_point(size=1.3) + 
    geom_linerange(size=1.1) + 
    geom_hline(yintercept = 0.5,alpha=0.3,col="red") +
    facet_wrap(~ Cntry) + 
    theme_bw() +
    ylim(0.3,0.7) +
    scale_x_continuous(breaks=xaxis_ticks,
                       labels = xaxis_lab) +
    geom_vline(xintercept=vertival_lines, col="grey", linetype="dashed", size=1) +
    coord_flip() +  xlab("Patients' Attributes") + ylab("Marginal Mean of Prioritization for Intensive Medical Care") + 
    scale_color_grey() +
    theme(legend.position = "none", text = element_text(size=14),
          strip.background =element_rect(fill="white"),
          #panel.grid.major.y =  element_blank(),
          panel.grid.minor.y =  element_blank())
    
  
  ggsave("maintext_figure1.pdf", height=9, width=12)
  
# Model AMCIE Affective Polarization ==========
  
  # Model
  est_lm_cond <- function(d,cuts_quantiles=F,p_cuts=c(0,0.33,0.66,1)){
    
    # Define cuts for rating
    if(cuts_quantiles) {
      cuts_rat <- quantile(d$resp_pview1 - d$resp_pview2,probs = p_cuts, na.rm = T)} 
    else {
      cuts_rat <- c(-10,-5,5,10)}
                   
    
    # Creat Diff
    d <- mutate(d, 
                drat = resp_pview1 - resp_pview2,
                dratc = cut(drat,cuts_rat))
    
    # Model Estimates
    m <- lm(y ~ partisanship1*dratc + partisanship2*dratc +
              age + arrival + chance + child + gender
            , data=d, weight=weight)
    clse_models <-  vcovCL(m, ~ ResponseId)
    
    # Create Estimates
    p <- coef(m)
    se <- sqrt(diag(clse_models))
    
    sel_pid1 <- grep("partisanship1",names(p),value = T)
    sel_pid2 <- grep("partisanship2",names(p),value = T)
    
    df_plot <- data.frame(
      "est" = c(p[sel_pid1[1]],p[sel_pid1[1]] + p[sel_pid1[2]],p[sel_pid1[1]] + p[sel_pid1[3]],  # effect among rest
                p[sel_pid2[1]],p[sel_pid2[1]] + p[sel_pid2[2]],p[sel_pid2[1]] + p[sel_pid2[3]]  # Among partisians party 2
      ),
      "se" = c(se[sel_pid1[1]],
               sqrt(clse_models[sel_pid1[1],sel_pid1[1]] + clse_models[sel_pid1[2],sel_pid1[2]] + 2*clse_models[sel_pid1[1],sel_pid1[2]]),
               sqrt(clse_models[sel_pid1[1],sel_pid1[1]] + clse_models[sel_pid1[3],sel_pid1[3]] + 2*clse_models[sel_pid1[1],sel_pid1[3]]),
               se[sel_pid2[1]],
               sqrt(clse_models[sel_pid2[1],sel_pid2[1]] + clse_models[sel_pid2[2],sel_pid2[2]] + 2*clse_models[sel_pid2[1],sel_pid2[2]]),
               sqrt(clse_models[sel_pid2[1],sel_pid2[1]] + clse_models[sel_pid2[3],sel_pid2[3]] + 2*clse_models[sel_pid2[1],sel_pid2[3]])),
      "rating" = c("Left Respondent","Middle","Right Respondent"),
      "pid" = c(rep("Right Party",3), rep("Left Party",3)))
    
    return(df_plot)
  } 
  
  # Estimation countrywise
  df_res <- df %>% 
    group_by(Cntry) %>%
    do(est_lm_cond(.))
  
  
  # Create Frame for plot
  df_plot <- filter(df_res)   %>%
    mutate(pat_pid = gsub("Party","PID", pid),
           resp_pid = gsub("Respondent ","",rating)) %>%
    mutate(label = paste(resp_pid,pat_pid)) %>%
    mutate(x_var =  case_when(label == "Right Respondent Right PID" ~ 1.4,
                              label == "Right Respondent Left PID" ~ 2.6,
                              
                              label == "Middle Right PID" ~ 4.4,
                              label == "Middle Left PID" ~ 5.6,
                              
                              label == "Left Respondent Right PID" ~ 7.4,
                              label == "Left Respondent Left PID" ~ 8.6,
                              
    )) %>%
    mutate(Cntry = factor(Cntry, levels = c("Countries combined",
                                            "USA","Germany","Italy","Poland","Brazil"))) 
  
  
  xaxis_ticks <- c(unique(df_plot$x_var), 2,5,8)
  xaxis_lab <- c("","","","","","","Right\n Party\n \n","Middle", "Left\n Party\n \n") 
  vertival_lines <- c(3.5,6.5)
  
  
  df_plot$pat_pid_lab <- case_when(
    df_plot$pat_pid == "No PID" ~ "No PID",
    df_plot$Cntry == "Countries combined" & df_plot$pat_pid == "Left PID" ~ "Patient Left PID",
    df_plot$Cntry == "Countries combined" & df_plot$pat_pid == "Right PID" ~ "Patient Right PID",
    df_plot$Cntry == "USA" & df_plot$pat_pid == "Left PID" ~ "Democrat",
    df_plot$Cntry == "USA" & df_plot$pat_pid == "Right PID" ~ "Republican",
    df_plot$Cntry == "Italy" & df_plot$pat_pid == "Left PID" ~ "Italia Viva",
    df_plot$Cntry == "Italy" & df_plot$pat_pid == "Right PID" ~ "Lega",
    df_plot$Cntry == "Germany" & df_plot$pat_pid == "Left PID" ~ "Green",
    df_plot$Cntry == "Germany" & df_plot$pat_pid == "Right PID" ~ "AfD",
    df_plot$Cntry == "Brazil" & df_plot$pat_pid == "Left PID" ~ "PSL",
    df_plot$Cntry == "Brazil" & df_plot$pat_pid == "Right PID" ~ "PT",
    df_plot$Cntry == "Poland" & df_plot$pat_pid == "Left PID" ~ "PO",
    df_plot$Cntry == "Poland" & df_plot$pat_pid == "Right PID" ~ "PiS",
    TRUE ~ df_plot$pat_pid
  )
  
  
  
  # Plot
  ggplot(df_plot,
         aes(x=x_var, y=est,ymin=est-1.96*se,ymax=est+1.96*se,label=pat_pid_lab, col=pat_pid)) +
    geom_pointrange(position = position_dodge2(0.5)) +
    geom_text(nudge_x = -0.4) +
    facet_wrap(~Cntry) +
    ylim(-0.3,0.3) +  
    #  geom_text(colour="darkgray", aes(y=0.2, label=resp_pid),  
    #            position=position_dodge(width=-0.9), col=gray) +
    geom_hline(yintercept=0, col="red",  alpha=0.3,size=1) +
    coord_flip(clip = "off") +  xlab("Respondents' Affective Polarization in favor of\n (Difference in Party Rating) \n") + 
    ylab("\nAMCIE (compared to patient with no PID)") + 
    scale_color_manual(values = c("black", "black", "black")) +
    scale_x_continuous(breaks=xaxis_ticks,
                       labels = xaxis_lab,limits = c(0.5,9.5)) +
    geom_vline(xintercept=vertival_lines, col="grey", linetype="dashed", size=1) +
    theme_bw() +
    theme(legend.position = "none", 
          text = element_text(size=19),
          axis.text.x = element_text(color = "black"),
          strip.background =element_rect(fill="white"),
          axis.text.y = element_text(color = "black", angle = 90,hjust = 0.5, vjust = 1),
          #panel.grid.major.y =  element_blank(),
          panel.grid.minor.y =  element_blank()
    ) 
  
  ggsave("maintext_figure2.pdf",width = 12, height = 9)  
  
  # Conditional marginal means
  est_mm_cond <- function(d,cuts_quantiles=F,p_cuts=c(0,0.33,0.66,1)){
    
    # Define cuts for rating
    if(cuts_quantiles) {
      cuts_rat <- quantile(d$resp_pview1 - d$resp_pview2,probs = p_cuts, na.rm = T)} 
    else {
      cuts_rat <- c(-10,-5,5,10)}
    
    # Creat Diff
    d <- mutate(d, 
                drat = resp_pview1 - resp_pview2,
                resp_pid = cut(drat, breaks = cuts_rat, labels = c("left","mid","right"))
                #resp_pid = cut(drat, c(-10,-5,5,10), labels = c("left","mid","right"))
    )
    
    # Model Estimates
    m <- lm(y ~ partisanship1*resp_pid + partisanship2*resp_pid +
              age + arrival + chance + child + gender
            , data=d, weight=weight)
    clse_models <-  vcovCL(m, ~ ResponseId)
    
    
    
    # Caclualte Marginal Means
    vars <- c("partisanship1","partisanship2", "resp_pid")
    df_res <- emmeans(m,vars, vcov.fun=clse_models)
    
    # Return
    return(as.data.frame(df_res))
  } 
  
  
  # Estimation countrywise
  df_res <- df %>% 
    group_by(Cntry) %>%
    do(est_mm_cond(.)) 
  
  # Create Frame for plot
  df_plot <- filter(df_res, 
                    #resp_pid!="mid",
                    !((partisanship1==1) & (partisanship2==1)))   %>%
    mutate(pat_pid = case_when(partisanship1 ==1 ~ "Right PID",
                               partisanship2 ==1 ~ "Left PID",
                               TRUE ~ "No PID")) %>%
    mutate(label = paste(resp_pid, pat_pid)) %>%
    mutate(x_var =  case_when(label == "right Right PID" ~ 1.2,
                              label == "right No PID" ~ 2,
                              label == "right Left PID" ~ 2.8,
                              
                              label == "mid Right PID" ~ 4.2,
                              label == "mid No PID" ~ 5,
                              label == "mid Left PID" ~ 5.8,
                              
                              label == "left Right PID" ~ 7.2,
                              label == "left No PID" ~ 8,
                              label == "left Left PID" ~ 8.8,
                              
    )) %>%
    mutate(Cntry = factor(Cntry, levels = c("Countries combined",
                                            "USA","Germany","Italy","Poland","Brazil"))) 
  
  
  xaxis_ticks <- unique(df_plot$x_var)
  xaxis_lab <- c("Left\n Party\n \n","","",
                 "Middle\n \n","","",
                 "Right\n Party\n \n","","") 
  vertival_lines <- c(3.5,6.5)
  
  
  df_plot$pat_pid_lab <- case_when(
    df_plot$pat_pid == "No PID" ~ "No PID",
    df_plot$Cntry == "Countries combined" & df_plot$pat_pid == "Left PID" ~ "Patient Left PID",
    df_plot$Cntry == "Countries combined" & df_plot$pat_pid == "Right PID" ~ "Patient Right PID",
    df_plot$Cntry == "USA" & df_plot$pat_pid == "Left PID" ~ "Democrat",
    df_plot$Cntry == "USA" & df_plot$pat_pid == "Right PID" ~ "Republican",
    df_plot$Cntry == "Italy" & df_plot$pat_pid == "Left PID" ~ "Italia Viva",
    df_plot$Cntry == "Italy" & df_plot$pat_pid == "Right PID" ~ "Lega",
    df_plot$Cntry == "Germany" & df_plot$pat_pid == "Left PID" ~ "Green",
    df_plot$Cntry == "Germany" & df_plot$pat_pid == "Right PID" ~ "AfD",
    df_plot$Cntry == "Brazil" & df_plot$pat_pid == "Left PID" ~ "PSL",
    df_plot$Cntry == "Brazil" & df_plot$pat_pid == "Right PID" ~ "PT",
    df_plot$Cntry == "Poland" & df_plot$pat_pid == "Left PID" ~ "PO",
    df_plot$Cntry == "Poland" & df_plot$pat_pid == "Right PID" ~ "PiS",
    TRUE ~ df_plot$pat_pid
  )
  
  
  df_plot$resp_pid <- case_when(df_plot$resp_pid == "right" ~ "Right\n Respondent",
                                df_plot$resp_pid == "left" ~ "Left Respondent",
                                df_plot$resp_pid == "mid" ~ "Middle\n Respondent",
  )
  
  
  
  # Plot
  ggplot(df_plot,
         aes(x=x_var, y=emmean,ymin=lower.CL,ymax=upper.CL,label=pat_pid_lab, col=pat_pid)) +
    geom_pointrange(position = position_dodge2(0.5)) +
    geom_text(nudge_x = -0.4) +
    facet_wrap(~Cntry) +
    ylim(0.3,0.7) +  
    #  geom_text(colour="darkgray", aes(y=0.2, label=resp_pid),  
    #            position=position_dodge(width=-0.9), col=gray) +
    geom_hline(yintercept=0.5, col="red",  alpha=0.1,size=1) +
    coord_flip(clip = "off") +  xlab("Respondents' Affective Polarization in favor of\n (Difference in Party Rating) \n") + 
    ylab("\nMarginal Mean of Prioritization for Intensive Medical Care") + 
    scale_color_manual(values = c("black", "black", "black")) +
    scale_x_continuous(breaks=xaxis_ticks,
                       labels = xaxis_lab,limits = c(0.5,9.5)) +
    geom_vline(xintercept=vertival_lines, col="grey", linetype="dashed", size=1) +
    theme_bw() +
    theme(legend.position = "none", 
          text = element_text(size=19),
          axis.text.x = element_text(color = "black"),
          strip.background =element_rect(fill="white"),
          axis.text.y = element_text(color = "black", angle = 90,hjust = 0.5, vjust = 1),
          #panel.grid.major.y =  element_blank(),
          panel.grid.minor.y =  element_blank()
    ) 
  
  ggsave("appendix_figure12.pdf",width = 12, height = 9)
  

  
# Model AMCIE Affective Polarization (quantiles) ==========
  
  # Estimation countrywise
  df_res <- df %>% 
    group_by(Cntry) %>%
    do(est_lm_cond(.,cuts_quantiles = T))
  
  
  # Create Frame for plot
  df_plot <- filter(df_res)   %>%
    mutate(pat_pid = gsub("Party","PID", pid),
           resp_pid = gsub("Respondent ","",rating)) %>%
    mutate(label = paste(resp_pid,pat_pid)) %>%
    mutate(x_var =  case_when(label == "Right Respondent Right PID" ~ 1.4,
                              label == "Right Respondent Left PID" ~ 2.6,
                              
                              label == "Middle Right PID" ~ 4.4,
                              label == "Middle Left PID" ~ 5.6,
                              
                              label == "Left Respondent Right PID" ~ 7.4,
                              label == "Left Respondent Left PID" ~ 8.6,
                              
    )) %>%
    mutate(Cntry = factor(Cntry, levels = c("Countries combined",
                                            "USA","Germany","Italy","Poland","Brazil"))) 
  
  
  xaxis_ticks <- c(unique(df_plot$x_var), 2,5,8)
  xaxis_lab <- c("","","","","","","Right\n Party\n \n","Middle", "Left\n Party\n \n") 
  vertival_lines <- c(3.5,6.5)
  
  
  df_plot$pat_pid_lab <- case_when(
    df_plot$pat_pid == "No PID" ~ "No PID",
    df_plot$Cntry == "Countries combined" & df_plot$pat_pid == "Left PID" ~ "Patient Left PID",
    df_plot$Cntry == "Countries combined" & df_plot$pat_pid == "Right PID" ~ "Patient Right PID",
    df_plot$Cntry == "USA" & df_plot$pat_pid == "Left PID" ~ "Democrat",
    df_plot$Cntry == "USA" & df_plot$pat_pid == "Right PID" ~ "Republican",
    df_plot$Cntry == "Italy" & df_plot$pat_pid == "Left PID" ~ "Italia Viva",
    df_plot$Cntry == "Italy" & df_plot$pat_pid == "Right PID" ~ "Lega",
    df_plot$Cntry == "Germany" & df_plot$pat_pid == "Left PID" ~ "Green",
    df_plot$Cntry == "Germany" & df_plot$pat_pid == "Right PID" ~ "AfD",
    df_plot$Cntry == "Brazil" & df_plot$pat_pid == "Left PID" ~ "PSL",
    df_plot$Cntry == "Brazil" & df_plot$pat_pid == "Right PID" ~ "PT",
    df_plot$Cntry == "Poland" & df_plot$pat_pid == "Left PID" ~ "PO",
    df_plot$Cntry == "Poland" & df_plot$pat_pid == "Right PID" ~ "PiS",
    TRUE ~ df_plot$pat_pid
  )
  
  
  
  # Plot
  ggplot(df_plot,
         aes(x=x_var, y=est,ymin=est-1.96*se,ymax=est+1.96*se,label=pat_pid_lab, col=pat_pid)) +
    geom_pointrange(position = position_dodge2(0.5)) +
    geom_text(nudge_x = -0.4) +
    facet_wrap(~Cntry) +
    ylim(-0.3,0.3) +  
    #  geom_text(colour="darkgray", aes(y=0.2, label=resp_pid),  
    #            position=position_dodge(width=-0.9), col=gray) +
    geom_hline(yintercept=0, col="red",  alpha=0.3,size=1) +
    coord_flip(clip = "off") +  xlab("Respondents' Affective Polarization in favor of\n (Difference in Party Rating) \n") + 
    ylab("\nAMCIE (compared to patient with no PID)") + 
    scale_color_manual(values = c("black", "black", "black")) +
    scale_x_continuous(breaks=xaxis_ticks,
                       labels = xaxis_lab,limits = c(0.5,9.5)) +
    geom_vline(xintercept=vertival_lines, col="grey", linetype="dashed", size=1) +
    theme_bw() +
    theme(legend.position = "none", 
          text = element_text(size=19),
          axis.text.x = element_text(color = "black"),
          strip.background =element_rect(fill="white"),
          axis.text.y = element_text(color = "black", angle = 90,hjust = 0.5, vjust = 1),
          #panel.grid.major.y =  element_blank(),
          panel.grid.minor.y =  element_blank()
    ) 
  
  ggsave("appendix_figure13.pdf",width = 12, height = 9)  
  
  
# Model AMCIE PID ==========
  
  # Model
  est_lm_cond <- function(d){
    
    # Model Estimates
    m <- lm(y ~ partisanship1*resp_pid1 + partisanship2*resp_pid2 +
              partisanship1*resp_pid2 + partisanship2*resp_pid1  +
              age + arrival + chance + child + gender
              , data=d, weight=weight)
    clse_models <-  vcovCL(m, ~ ResponseId)

    # Create Estimates
    p <- coef(m)
    se <- sqrt(diag(clse_models))
    
    df_plot <- data.frame(
      "est" = c(p["partisanship1"], p["partisanship2"], # effect among rest
                p["partisanship1"] + p["partisanship1:resp_pid1"], p["partisanship2"] + p["resp_pid1:partisanship2"], # Among partisians party 1
                p["partisanship1"] + p["partisanship1:resp_pid2"], p["partisanship2"] + p["partisanship2:resp_pid2"]  # Among partisians party 2
                ),
      "se" = c(se["partisanship1"],se["partisanship2"],
               sqrt(clse_models["partisanship1","partisanship1"] + clse_models["partisanship1:resp_pid1","partisanship1:resp_pid1"] + 2*clse_models["partisanship1","partisanship1:resp_pid1"]),
               sqrt(clse_models["partisanship2","partisanship2"] + clse_models["resp_pid1:partisanship2","resp_pid1:partisanship2"] + 2*clse_models["partisanship1","resp_pid1:partisanship2"]),
               sqrt(clse_models["partisanship1","partisanship1"] + clse_models["partisanship1:resp_pid2","partisanship1:resp_pid2"] + 2*clse_models["partisanship1","partisanship1:resp_pid2"]),
               sqrt(clse_models["partisanship2","partisanship2"] + clse_models["partisanship2:resp_pid2","partisanship2:resp_pid2"] + 2*clse_models["partisanship1","partisanship2:resp_pid2"])),
      "party" = c("Patient Right Partisan ","Left Partisan Patient"),
      "pid" = c(rep("None/Other Party",2),rep("Right Party",2), rep("Left Party",2)))
    
    return(df_plot)
  } 

  # Estimation countrywise
  df_res <- df %>% 
    group_by(Cntry) %>%
    do(est_lm_cond(.))
  
  # Prepare Df
  # Create Frame for plot
  df_plot <- filter(df_res)   %>%
    mutate(pat_pid = gsub("Partisan | Patient ","", party),
           resp_pid = gsub("","",pid)) %>%
    mutate(label = paste(resp_pid,pat_pid)) %>%
    mutate(x_var =  case_when(label == "Right Party Patient Right " ~ 1.2,
                              label == "Right Party Left Patient" ~ 1.8,
                              
                              label == "None/Other Party Patient Right " ~ 3.2,
                              label == "None/Other Party Left Patient" ~ 3.8,
                              
                              label == "Left Party Patient Right " ~ 5.2,
                              label == "Left Party Left Patient" ~ 5.8
                              
    )) %>%
    mutate(Cntry = factor(Cntry, levels = c("Countries combined",
                                            "USA","Germany","Italy","Poland","Brazil"))) 
  
  
  xaxis_ticks <- c(unique(df_plot$x_var),1.5,3.5,5.5)
  xaxis_lab <- c("","","","","","","Right\n Party\n \n","None/Other\n Party","Left\n Party\n \n") 
  vertival_lines <- c(2.5,4.5)
  
  
  df_plot$pat_pid_lab <- case_when(
    df_plot$pat_pid == "No PID" ~ "No PID",
    df_plot$Cntry == "Countries combined" & df_plot$pat_pid == "Left Patient" ~ "Patient Left PID",
    df_plot$Cntry == "Countries combined" & df_plot$pat_pid == "Patient Right " ~ "Patient Right PID",
    df_plot$Cntry == "USA" & df_plot$pat_pid == "Left Patient" ~ "Democrat",
    df_plot$Cntry == "USA" & df_plot$pat_pid == "Patient Right " ~ "Republican",
    df_plot$Cntry == "Italy" & df_plot$pat_pid == "Left Patient" ~ "Italia Viva",
    df_plot$Cntry == "Italy" & df_plot$pat_pid == "Patient Right " ~ "Lega",
    df_plot$Cntry == "Germany" & df_plot$pat_pid == "Left Patient" ~ "Green",
    df_plot$Cntry == "Germany" & df_plot$pat_pid == "Patient Right " ~ "AfD",
    df_plot$Cntry == "Brazil" & df_plot$pat_pid == "Left Patient" ~ "PSL",
    df_plot$Cntry == "Brazil" & df_plot$pat_pid == "Patient Right " ~ "PT",
    df_plot$Cntry == "Poland" & df_plot$pat_pid == "Left Patient" ~ "PO",
    df_plot$Cntry == "Poland" & df_plot$pat_pid == "Patient Right " ~ "PiS",
    TRUE ~ df_plot$pat_pid
  )
  
  
  
  # Plot
  ggplot(df_plot,
         aes(x=x_var, y=est,ymin=est-1.96*se,ymax=est+1.96*se,label=pat_pid_lab, col=pat_pid)) +
    geom_pointrange(position = position_dodge2(0.5)) +
    geom_text(nudge_x = -0.2) +
    facet_wrap(~Cntry) +
    ylim(-0.5,0.5) +  
    #  geom_text(colour="darkgray", aes(y=0.2, label=resp_pid),  
    #            position=position_dodge(width=-0.9), col=gray) +
    geom_hline(yintercept=0, col="red",  alpha=0.3,size=1) +
    coord_flip(clip = "off") +  xlab("Partisan Identification with \n") + 
    ylab("\nAMCIE (compared to patient with no PID)") + 
    scale_color_manual(values = c("black", "black", "black")) +
    scale_x_continuous(breaks=xaxis_ticks,
                       labels = xaxis_lab,limits = c(0.5,6.5)) +
    geom_vline(xintercept=vertival_lines, col="grey", linetype="dashed", size=1) +
    theme_bw() +
    theme(legend.position = "none", 
          text = element_text(size=19),
          axis.text.x = element_text(color = "black"),
          strip.background =element_rect(fill="white"),
          axis.text.y = element_text(color = "black", angle = 90,hjust = 0.5, vjust = 1),
          #panel.grid.major.y =  element_blank(),
          panel.grid.minor.y =  element_blank()
    ) 
  
  ggsave("appendix_figure16.pdf",width = 12, height = 9) 
 

  
  
# Model AMCIE Affective Polarization (Linear) ==========
  
  # Model
  est_lm_cond <- function(d){
    
    # Creat Diff
    d <- mutate(d, 
                drat = resp_pview1 - resp_pview2)
    
    # Model Estimates
    m <- lm(y ~ partisanship1*drat + partisanship2*drat +
              age + arrival + chance + child + gender
            , data=d, weight=weight)
    clse_models <-  vcovCL(m, ~ ResponseId)
    
    # Create Estimates
    p <- coef(m)
    se <- sqrt(diag(clse_models))
    
    sel_pid1 <- grep("partisanship1",names(p),value = T)
    sel_pid2 <- grep("partisanship2",names(p),value = T)
    
    df_plot <- data.frame(
      "est" = c(p[sel_pid1[1]],p[sel_pid1[2]],  # effect among rest
                p[sel_pid2[1]],p[sel_pid2[2]],p["drat"]   # Among partisians party 2
      ),
      "se" = c(c(se[sel_pid1[1]],se[sel_pid1[2]],  # effect among rest
                 se[sel_pid2[1]],se[sel_pid2[2]],se["drat"]   # Among partisians party 2
      )),
      "par" = c("Right Partisanship","Right Partisanship X Affective Polarization",
                "Left Partisanship","Left Partisanship X Affective Polarization","Affective Polarization"))
    
    return(df_plot)
  } 
  
  # Estimation countrywise
  df_res <- df %>% 
    group_by(Cntry) %>%
    do(est_lm_cond(.))
  
  
  # Prepare Df
  df_plot <- df_res %>%
    mutate(Cntry = factor(Cntry, levels = c("Countries combined","USA","Germany","Italy","Poland","Brazil")))
  
  
  # Figure
  ggplot(df_plot) + 
    geom_hline(yintercept = 0,alpha=0.3,col="red") +
    geom_pointrange(aes(x=par,y=est,ymin=est-2*se,ymax=est+2*se),size=.8) +
    theme_bw() +
    coord_flip() + facet_wrap(~ Cntry) +
    scale_color_grey(name="Rating Differential") +  theme(text = element_text(size=18)) +
    xlab("") + ylab("ACME") +
    theme(legend.position = "bottom", 
          text = element_text(size=12),
          #panel.grid.major.y =  element_blank(),
          panel.grid.minor.y =  element_blank())
  
  
  ggsave("appendix_figure17.pdf",width = 12, height = 9)  
  
  
# Model AMCIE Sociodemgraphics ==========
  
  # Model
  est_lm_cond_soc <- function(d, var_resp="male", var_par="gender"){
    
    # Formula
    if(var_par == "age_patient"){
      form <- as.formula(paste("y ~ partisanship1 + partisanship2 +
              age_patient + arrival + chance + child + gender + ", var_resp, "*" ,var_par))  
    } else {
      form <- as.formula(paste("y ~ partisanship1 + partisanship2 +
              age + arrival + chance + child + gender + ", var_resp, "*" ,var_par))  
    }
  
    # Model Estimates
    m <- lm(form, data=d, weight=weight)
    clse_models <-  vcovCL(m, ~ ResponseId)
    
    # Create Estimates
    p <- coef(m)
    sel_var <- grep(var_par,names(p),value = T)
    sel_int <- sel_var[grepl(":",sel_var)]
    sel_var <- sel_var[!grepl(":",sel_var) & sel_var != var_resp]
      
      
    # Get Coeffients
    d_eff <- p[sel_var]
    cmbd_eff <- p[sel_var] + p[sel_int]
    int_eff <- p[sel_int]
    
    # Standard Error
    d_se <- sqrt(diag(clse_models))[sel_var]
    int_se <- sqrt(diag(clse_models))[sel_int]
    cmbd_se <- sqrt(clse_models[sel_var,sel_var] + clse_models[sel_int,sel_int] + 2* clse_models[sel_var,sel_int])
    
    # Data.frame
    df_plot <- data.frame(
      "effect" = c("cond_eff","cond_eff","inter_eff"),
      "resp_attr" = c(0,1,0),
      "eff" = c(d_eff,cmbd_eff,int_eff),
      "se" =  c(d_se,cmbd_se,int_se),
      "val" = c(sel_var, paste(sel_var,"+",sel_int), sel_int)) %>%
      mutate("p" = round(1-pnorm(abs(eff/se)),2))
    
    return(df_plot)
  } 
  
  
  # Estimation Gender countrywise
  df_res_gender <- df %>% 
    group_by(Cntry) %>%
    do(est_lm_cond_soc(.)) 
  
  # Estimation Children countrywise
  df_res_child <- df %>% 
    group_by(Cntry) %>%
    do(est_lm_cond_soc(.,var_resp = "children_household",var_par = "child")) 
  
 
  # Estimation Age countrywise
  df_res_age <- df %>% 
    mutate(age_resp = ifelse(age_resp > 60, 1,0),
           age_patient = ifelse(age > 60,1,0)) %>% 
    group_by(Cntry) %>%
    do(est_lm_cond_soc(.,var_resp = "age_resp",var_par = "age_patient")) 
  

  # Combine Estimates  
  df_res_socio <- bind_rows(df_res_gender,df_res_child,df_res_age)

  filter(df_res_socio, val == "genderMale + genderMale:male")
  

  # Prepare Df
  df_plot <- filter(df_res_socio, effect == "inter_eff") %>%
    mutate(Cntry = factor(Cntry, levels = c("Countries combined","USA","Germany","Italy","Poland","Brazil")),
           val = as.factor(case_when(val == "genderMale:male" ~ "Male Patient X Male Respond.",
                           val == "childYes:children_household" ~ "Patient w Children X Resp. w Children",
                           val == "age_patient:age_resp" ~ "Patient above Age 60 X Resp. above Age 60")
           ))
  
  # Figure
  ggplot(df_plot) + 
    geom_hline(yintercept = 0,alpha=0.3,col="red") +
    geom_pointrange(aes(x=reorder(val,eff),y=eff,ymin=eff-1.67*se,ymax=eff+1.67*se),
                    position = position_dodge(0.4), size=.8) +
    theme_bw() +
    coord_flip() + facet_wrap(~ Cntry) +
    scale_color_grey(name="Respondent PID") +  theme(text = element_text(size=18)) +
    xlab("") + ylab("ACME") +
    theme(legend.position = "none", text = element_text(size=12),
          #panel.grid.major.y =  element_blank(),
          panel.grid.minor.y =  element_blank())
  
  # Save
  ggsave("appendix_figure15.pdf",width = 12, height = 9)  
  
  
# Test Random Forest Variable Importance ================
  

  # Estimation Function AMCE
  est_randomforest <- function(df){
    df_model <- df %>%
      select(y,age , arrival , chance , child , job , gender , partisanship1 , partisanship2) %>%
      na.omit()
    
    X <- model.matrix(y ~ age + arrival + chance + child + job + gender + partisanship1 + partisanship2,df_model)
    y <- factor(df_model$y)
    fit <- randomForest(X[,-1],y)
    var_import <- importance(fit)
    
    df_var_import <- data.frame(var_import)
    df_var_import$var <- rownames(df_var_import)
    rownames(df_var_import) <- NULL
    
    return(df_var_import)
  } 

  # Results
  df_names <- data.frame(
    var = c("age42","age61", "age76",
            "jobNurse","jobPhysician","jobProfessor","jobUnemployed",                
            "arrivalSecond","arrivalSame time",
            "chance50%","chance80%",
            "partisanship1","partisanship2",
            "genderMale",
            "childYes"),
    attr   = c("Age","Age","Age", 
               "Job","Job","Job","Job",                  
               "Arrival","Arrival",
               "Chance of Survival","Chance of Survival", 
               "Partisanship","Partisanship",
               "Gender",
               "Children"),
    levl = c("42","61", "76",
             "Nurse","Physician","Professor","Unemployed",               
             "Second","Same time",
             "50%","80%",
             "Right","Left",
             "Male",
             "Yes"))
  
  # Estimation countrywise
  df_res <- df %>% 
    group_by(Cntry) %>%
    do(est_randomforest(.))

  
  df_plot <- left_join(df_res,df_names) %>%
    mutate(Cntry = factor(Cntry, levels = c("Countries combined","USA","Germany","Italy","Poland","Brazil")),
           var_lab <- paste(attr, levl, sep=" "),
           var_lab = reorder_within(var_lab,MeanDecreaseGini,Cntry, sep= " "))
  

  # Figure
  ggplot(df_plot) + 
    geom_point(aes(x=var_lab,y=MeanDecreaseGini), size=.8) +
    theme_bw() +
    coord_flip() + facet_wrap(~ Cntry, scales = "free") +
    scale_color_grey(name="Respondent PID") +  theme(text = element_text(size=18)) +
    xlab("") + ylab("Mean Decrease Gini") +
    theme(legend.position = "none", text = element_text(size=12),
          #panel.grid.major.y =  element_blank(),
          panel.grid.minor.y =  element_blank())

  
  # Save
  ggsave("appendix_figure21.pdf",width = 12, height = 9)  
  
# Estimation among Respondents who pass Attention Check ================

  # Estimation countrywise
  df_res <- df %>% 
    filter(cjattcheck) %>%
    group_by(Cntry) %>%
    do(est_lm(.))
  
  df %>% 
    filter(cjattcheck) %>%
    group_by(Cntry) %>%
    summarise(n())

  # List of ACMEs
  df_names <- data.frame(
    names = c("age21","age42","age61", "age76",
              "jobNurse","jobPhysician","jobProfessor","jobUnemployed","jobCook",                
              "arrivalFirst","arrivalSecond","arrivalSame time",
              "chance20%","chance50%","chance80%",
              "partisanship1","partisanship2","partisanship3",
              "genderFemale","genderMale",
              "childYes","childNo"),
    attr   = c("Age","Age","Age", "Age",
               "Job","Job","Job","Job","Job",                  
               "Arrival","Arrival","Arrival",
               "Chance of Survival","Chance of Survival", "Chance of Survival",
               "Partisanship","Partisanship","Partisanship",
               "Gender","Gender",
               "Children","Children"),
    levl = c("21","42","61", "76",
             "Nurse","Physician","Professor","Unemployed","Cook",                  
             "First","Second","Same time",
             "20%","50%","80%",
             "Right","Left","None",
             "Female","Male",
             "Yes","No"))
  
  df_names <- crossing(df_names, data.frame("Cntry"=unique(df_res$Cntry)))
  
  # Create Baseline
  df_plot <- left_join(df_names, df_res) %>%
    mutate(est = ifelse(is.na(est),0,est),
           se = ifelse(is.na(se),0,se),
           baseline = ifelse(est == 0 & se ==0,"Baseline","ACME"),
           levl = as.factor(levl), 
           Cntry = factor(Cntry, levels = c("Countries combined","USA","Germany","Italy","Poland","Brazil"))) 
  
  df_plot$label <- as.factor(interaction(df_plot$levl, df_plot$attr))
  df_plot$position <- case_when(df_plot$label == "21.Age" ~ 1,
                                df_plot$label == "42.Age" ~ 2,
                                df_plot$label == "61.Age" ~ 3,
                                df_plot$label == "76.Age" ~ 4,
                                
                                df_plot$label == "First.Arrival" ~ 6,
                                df_plot$label == "Same time.Arrival" ~ 7,
                                df_plot$label == "Second.Arrival" ~ 8,
                                
                                df_plot$label == "20%.Chance of Survival" ~ 10,
                                df_plot$label == "50%.Chance of Survival" ~ 11,
                                df_plot$label == "80%.Chance of Survival" ~ 12,
                                
                                df_plot$label == "No.Children" ~ 14,
                                df_plot$label == "Yes.Children" ~ 15,
                                
                                df_plot$label == "Female.Gender" ~ 17,
                                df_plot$label == "Male.Gender" ~ 18,
                                
                                df_plot$label == "Cook.Job" ~ 20,
                                df_plot$label == "Nurse.Job" ~ 21,
                                df_plot$label == "Physician.Job" ~ 22,
                                df_plot$label == "Professor.Job" ~ 23,
                                df_plot$label == "Unemployed.Job" ~ 24,
                                
                                df_plot$label == "None.Partisanship" ~ 26,
                                df_plot$label == "Right.Partisanship" ~ 27,
                                df_plot$label == "Left.Partisanship" ~ 28)
  
  
  xaxis_ticks <- unique(df_plot$position)
  xaxis_lab <- gsub("\\."," ", unique(df_plot$label))
  
  vertival_lines <- c(5,9,13,16,19,25)
  
  
  # Figure
  ggplot(df_plot,aes(x=position,y=est,ymin=est-1.96*se,ymax=est+1.96*se, 
                     col=baseline)) + 
    geom_hline(yintercept = 0,alpha=0.3,col="red") +
    facet_wrap(~ Cntry) + theme_bw() +
    geom_pointrange(size=0.7) +
    scale_x_continuous(breaks=xaxis_ticks,
                       labels = xaxis_lab) +
    geom_vline(xintercept=vertival_lines, col="grey", linetype="dashed", size=1) +
    coord_flip() +  xlab("") + ylab("ACME") + 
    scale_color_grey() +
    theme(legend.position = "none", text = element_text(size=12),
          #panel.grid.major.y =  element_blank(),
          panel.grid.minor.y =  element_blank())
  
  ggsave("appendix_figure18.pdf", height=9, width=12) 
  
  
  # 